

### This script reproduces the main analyses. ###

library(sjPlot)
library(sjlabelled)
library(tidyverse)
library(data.table)
library(glm.predict)

# set working directory to folder with all scripts and data
setwd("/Users/janmattidollbaum/Dropbox/01_Academic_Texts/2020_Opposition_Personality/CPS/Code_and_Data")

##### load data #####

data1 <- fread("data1.csv")
data2 <- fread("data2.csv")
pooled <- fread("pooled.csv")
BYdata <- fread("BYdata.csv")

##### Hypothesis 1 #####

### Table G1
mH1_set1_1 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST,
              family = "binomial", data = data1)

mH1_set1_2 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                age + edu + female + wealth + moscow,
              family = "binomial", data = data1)

mH1_set1_3 <- glm(polact1 ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow + putin_sup,
               family = "binomial", data = data1)

H1_models <- list(mH1_set1_1, mH1_set1_2, mH1_set1_3)
H1_modelnames <- c("1", "2", "3")

tab_model(H1_models,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_G1.html")

### marginal effects

predictH1a_female <- ggeffects::ggeffect(mH1_set1_2, terms = "female [0,1]")
predictH1a_female$variable <- "female"
predictH1a_extra <- ggeffects::ggeffect(mH1_set1_2, terms = "extraST [meansd]")
predictH1a_extra$variable <- "extraversion"
predictH1a_moscow <- ggeffects::ggeffect(mH1_set1_2, terms = "moscow")
predictH1a_moscow$variable <- "Moscow"

### Figure 1

pH1_extra <- ggplot(predictH1a_extra[-c(2),]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of activism") +
  ylim(0.3, 0.7) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH1_extra

pH1_female <- ggplot(predictH1a_female) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("") +
  ylim(0.3, 0.7) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH1_female

pH1_moscow <- ggplot(predictH1a_moscow) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("") +
  ylim(0.3, 0.7) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH1_moscow

ggsave("Figure_1.pdf", gridExtra::arrangeGrob(pH1_extra, pH1_female, pH1_moscow, ncol = 3),
       width = 7, height = 3)


##### Hypothesis 2 #####

### Table G2

mH2_set1_1 <- glm(navact ~ openST + conscST + extraST + agreeST + neurST,
              family = "binomial", data = data2)

mH2_set1_2 <- glm(navact ~ openST + conscST + extraST + agreeST + neurST +
                age + edu + female + wealth + moscow,
              family = "binomial", data = data2)

mH2_set1_3 <- MASS::polr(factor(navact_scale) ~ openST + conscST + extraST + agreeST + neurST,
               method = "logistic", data = data2)

mH2_set1_4 <- MASS::polr(factor(navact_scale) ~ openST + conscST + extraST + agreeST + neurST +
                 age + edu + female + wealth + moscow,
               method = "logistic", data = data2)

H2_models <- list(mH2_set1_1, mH2_set1_2, mH2_set1_3, mH2_set1_4)
H2_modelnames <- c("1", "2", "3", "4")

tab_model(H2_models,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_G2.html")

### Figure 2

## marginal effects

predictH2_female <- ggeffects::ggeffect(mH2_set1_2, terms = "female [0,1]")
predictH2_female$variable <- "female"
predictH2_extra <- ggeffects::ggeffect(mH2_set1_2, terms = "extraST [meansd]")
predictH2_extra$variable <- "extraversion"
predictH2_moscow <- ggeffects::ggeffect(mH2_set1_2, terms = "moscow")
predictH2_moscow$variable <- "Moscow"

# plot
pH2_extra <- ggplot(predictH2_extra[-c(2),]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of offline activism") +
  ylim(0.1, 0.45) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH2_extra

pH2_female <- ggplot(predictH2_female) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("") +
  ylim(0.1, 0.45) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH2_female

pH2_moscow <- ggplot(predictH2_moscow) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("") +
  ylim(0.1, 0.45) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH2_moscow

ggsave("Figure_2.pdf", gridExtra::arrangeGrob(pH2_extra, pH2_female, pH2_moscow, ncol = 3),
       width = 7, height = 3)


### Figure 3

predictH2a <- data.frame(basepredict.polr(mH2_set1_4, c(0,0,-2,0,0,
                                                    mean(data2$age, na.rm = T),
                                                    mean(data2$edu, na.rm = T),
                                                    mean(data2$female, na.rm = T),
                                                    mean(data2$wealth, na.rm = T),
                                                    mean(data2$moscow, na.rm = T))))
predictH2a$x <- "-2sd"

predictH2b <- data.frame(basepredict.polr(mH2_set1_4, c(0,0,-1,0,0,
                                                    mean(data2$age, na.rm = T),
                                                    mean(data2$edu, na.rm = T),
                                                    mean(data2$female, na.rm = T),
                                                    mean(data2$wealth, na.rm = T),
                                                    mean(data2$moscow, na.rm = T))))
predictH2b$x <- "-1sd"

predictH2c <- data.frame(basepredict.polr(mH2_set1_4, c(0,0,0,0,0,
                                                    mean(data2$age, na.rm = T),
                                                    mean(data2$edu, na.rm = T),
                                                    mean(data2$female, na.rm = T),
                                                    mean(data2$wealth, na.rm = T),
                                                    mean(data2$moscow, na.rm = T))))
predictH2c$x <- "mean"

predictH2d <- data.frame(basepredict.polr(mH2_set1_4, c(0,0,1,0,0,
                                                    mean(data2$age, na.rm = T),
                                                    mean(data2$edu, na.rm = T),
                                                    mean(data2$female, na.rm = T),
                                                    mean(data2$wealth, na.rm = T),
                                                    mean(data2$moscow, na.rm = T))))
predictH2d$x <- "+1sd"

predictH2e <- data.frame(basepredict.polr(mH2_set1_4, c(0,0,2,0,0,
                                                    mean(data2$age, na.rm = T),
                                                    mean(data2$edu, na.rm = T),
                                                    mean(data2$female, na.rm = T),
                                                    mean(data2$wealth, na.rm = T),
                                                    mean(data2$moscow, na.rm = T))))
predictH2e$x <- "+2sd"

predictH2 <- rbind(predictH2a, predictH2b, predictH2c, predictH2d, predictH2e)
predictH2$outcome <- rep(c("none", "sporadic", "regular"), 5)
colnames(predictH2)[2:3] <- c("conf.min", "conf.max") 

predictH2$x <- factor(predictH2$x, levels = c("-2sd", "-1sd", "mean", "+1sd", "+2sd"))
predictH2$outcome <- factor(predictH2$outcome, levels = c("none", "sporadic", "regular"))

p.H2_polr <- ggplot(predictH2) +
  geom_pointrange(aes(x = factor(x), y = mean,
                      ymin = conf.min, ymax = conf.max)) +
  facet_wrap(factor(predictH2$outcome)) +
  xlab("extraversion") +
  ylab("predicted probability of campaign work intensity") +
  theme_bw(); p.H2_polr

ggsave("Figure_3.pdf", width = 9, height = 4.5)


##### Hypothesis 3 #####

# the variable 'subset' is constructed based on Table 2 in the main text, with
# subset==1 as the active oppositionists (i.e. the cases) and subset==0 the 
# non-active regime supporters (i.e. the controls)

### Table G3

# main models with 2013 data
mH3_set1_1 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST,
                 family = "binomial", data = data1)

mH3_set1_2 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                   age + edu + female + wealth + moscow,
                 family = "binomial", data = data1)

# main models with pooled data
mH3_set1_3 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST,
                 family = "binomial", data = pooled)

mH3_set1_4 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                   age + edu + female + wealth + moscow,
                 family = "binomial", data = pooled)

tab_model(mH3_set1_1, mH3_set1_2, mH3_set1_3, mH3_set1_4,
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_G3.html")

### Figure 4

# marginal effects
predictH3a_ex <- ggeffects::ggeffect(mH3_set1_2, terms = c("extraST [meansd]"))
predictH3a_ag <- ggeffects::ggeffect(mH3_set1_2, terms = c("agreeST [meansd]"))
predictH3a_female <- ggeffects::ggeffect(mH3_set1_2, terms = c("female [0,1]"))
predictH3a_moscow <- ggeffects::ggeffect(mH3_set1_2, terms = c("moscow [0,1]"))

predictH3b_ex <- ggeffects::ggeffect(mH3_set1_4, terms = c("extraST [meansd]"))
predictH3b_ag <- ggeffects::ggeffect(mH3_set1_4, terms = c("agreeST [meansd]"))
predictH3b_female <- ggeffects::ggeffect(mH3_set1_4, terms = c("female [0,1]"))
predictH3b_moscow <- ggeffects::ggeffect(mH3_set1_4, terms = c("moscow [0,1]"))


# plots for 2013 survey
pH3a_extra <- ggplot(predictH3a_ex[-c(2),]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_extra

pH3a_agree <- ggplot(predictH3a_ag[-c(2),]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("agreeableness") +
  ylab("\n") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_agree

pH3a_female <- ggplot(predictH3a_female) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("\n") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_female

pH3a_moscow <- ggplot(predictH3a_moscow) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("\n") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_moscow

ggsave("Figure_4a.pdf", gridExtra::grid.arrange(pH3a_extra, pH3a_agree, pH3a_female, pH3a_moscow, 
                               ncol = 4, top = "2013 survey"),
       width = 10, height = 3)


# plots for pooled data
pH3b_extra <- ggplot(predictH3b_ex[-c(2),]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  ylim(0, 0.85) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_extra

pH3b_agree <- ggplot(predictH3b_ag[-c(2),]) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("agreeableness") +
  ylab("\n") +
  ylim(0, 0.85) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_agree

pH3b_female <- ggplot(predictH3b_female) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("\n") +
  ylim(0, 0.85) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_female

pH3b_moscow <- ggplot(predictH3b_moscow) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("\n") +
  ylim(0, 0.85) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_moscow

ggsave("Figure_4b.pdf", gridExtra::grid.arrange(pH3b_extra, pH3b_agree, pH3b_female, pH3b_moscow, 
                               ncol = 4, top = "pooled 2013/Navalny"),
       width = 10, height = 3)


##### Hypothesis 4 #####

### Table K1

m4_set1_1 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                  age + edu + female + wealth + moscow +
                  extraST*agreeST,
                family = "binomial", data = data1)

m4_set1_2 <- glm(subset ~ openST + conscST + extraST + agreeST + neurST +
                  age + edu + female + wealth + moscow + 
                  extraST*agreeST,
                family = "binomial", data = pooled)

tab_model(m4_set1_1, m4_set1_2, 
          show.se = TRUE, show.ci = FALSE, 
          p.style = "stars", collapse.se = TRUE, show.ci50 = TRUE, show.aic = TRUE,
          file = "Table_K1.html")

### Figure 5

# prepare marginal effects tables
me.H4a_r <- ggeffects::ggeffect(m4_set1_1, terms = c("extraST", "agreeST"))
me.H4b_r <- ggeffects::ggeffect(m4_set1_2, terms = c("extraST", "agreeST"))

# merge marginal effects tables
me.H4a_r$group2 <- as.factor(rep.int(c("-1 sd", "mean", "+1 sd"), nrow(me.H4a_r)/3))
me.H4a_r$sample <- "2013 survey"

me.H4b_r$group2 <- as.factor(rep.int(c("-1 sd", "mean", "+1 sd"), nrow(me.H4b_r)/3))
me.H4b_r$sample <- "pooled: 2013/Navalny"

me.H4_r <- rbind(me.H4a_r, me.H4b_r)

# plot
me.H4_r$group2 <- factor(me.H4_r$group2, levels = c("-1 sd", "mean", "+1 sd"))
me.H4_r$sample <- factor(me.H4_r$sample, 
                         levels = c("2013 survey",
                                    "pooled: 2013/Navalny"))

ggplot(me.H4_r[me.H4_r$group2!="mean",]) +
  geom_smooth(aes(x=x, y=predicted, color = group2),
              size=.5, se = F) +
  geom_ribbon(aes(x=x, ymin = conf.low, ymax = conf.high, fill = group2), 
              alpha = .2) +
  scale_color_discrete(name="agreeableness") +
  scale_fill_discrete(name="agreeableness") +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  facet_wrap(~sample, 
             ncol = 3) +
  theme_bw()

ggsave("Figure_5.pdf", width = 7, height = 4)
